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We describe a simulation method for the accurate study of the equilibrium freezing properties of 
polydisperse fluids under the experimentally relevant condition of fixed polydispersity. The approach 
is based on the phase switch Monte Carlo method of Wilding and Bruce [Phys. Rev. Lett. 85, 
5138 (2000)]. This we have generalized to deal with particle size polydispersity by incorporating 
updates which alter the diameter a of a particle, under the control of a distribution of chemical 
potential differences p,(a). Within the resulting isobaric semi-grand canonical ensemble, we detail 
how to adapt p,(a) and the applied pressure such as to study coexistence, whilst ensuring that the 
ensemble averaged density distribution p(a) matches a fixed functional form. Results are presented 
for the effects of small degrees of polydispersity on the solid-liquid transition of soft spheres. 



PACS numbers: 

I. INTRODUCTION AND BACKGROUND 

Our ability to accurately determine the phase coex- 
istence properties of model systems via computer sim- 
ulation has increased greatly in recent years. Progress 
has come in the form of specialized Monte Carlo meth- 
ods that are tailor-made for tackling the core problem 
of coexistence, namely the comparison of the statistical 
weights of the disjoint regions of configuration space as- 
sociated with competing phases. To achieve this one gen- 
erally aims to engineer a sampling path which connects 
the phases in question, allowing each to be visited re- 
peatedly in the course of a single simulation run. The 
reward for doing so comes in the form of direct, accurate 
and transparent measurements of free energy differences 
and coexistence parameters 

In practice, however, the most obvious choice of an 
inter-phase sampling path, namely one routed via mixed- 
phase (ie. interfacial) configurations may not always be 
the most straightforward to negotiate. Typically such 
regions of configuration space are characterized by large 
free energy barriers, or by a 'pitted' free energy landscape 

- features which hinder efficient sampling. Contemporary 
strategies seek to either surmount or circumvent these 
impediments (for a review see [l|). One such scheme 

- which is in principle rather general - is phase switch 
Monte Carlo (PSMC). Its key feature is a phase space 
'leap' that directly maps a pure phase configuration 
of one phase onto a pure phase configuration of another. 
As such it completely avoids mixed-phase configurations 
and any attendant sampling difficulties. 

In the context of fluid-solid coexistence, PSMC per- 
mits efficient two-phase sampling and hence the direct de- 
termination of equilibrium freezing parameters and their 
uncertainties. This contrasts with simulation approaches 
which seek to emulate physical crystal nucleation pro- 
cesses, which commonly encounter significant problems, 
principally a large degree of metastability of the phases, 
extended timescales for crystallization and a tendency 
for the crystals formed to exhibit vacancies, interstitials 
and other defects. To date, PSMC has been successfully 
applied to the freezing of hard spheres, softly repulsive 



spheres, and the Lennard- Jones fluid [!, 0, H, H[. 

The present paper is concerned with extending the ap- 
plicability of PSMC to the study of the freezing prop- 
erties of model polydisperse fluids. Polydispersity is the 
feature, pervasive in soft condensed matter systems such 
as colloids and liquid crystals, whereby the constituent 
particles are not all identical but instead exhibit variation 
in terms of some attribute {a, say) such as the particle 
size or shape, etc. It leads to phase behaviour that is 
considerably richer in both variety and form than that of 
monodisperse systems 0- 

Typically one describes the polydispersity of a given 
system in terms of a density distribution p(°\a) which 
counts the number of particles per unit volume having 
a in the range a ... a + da. In most real polydisperse 
systems the form of p(°\a) is fixed by the synthesis of 
the fluid, and only its scale can change according to the 
degree of dilution of the system. Accordingly, one writes 
p(°)(cr) = n^f(a) where f(a) is a normalized fixed shape 
function and n<°> is the overall number density. Varying 
rv- ' corresponds to scanning a "dilution line" of the sys- 
tem. 

At coexistence, the particles of the various a values 
will be partitioned unequally between the phases. This is 
the phenomenon of "fractionation" which underpins the 
opulence of polydisperse phase behaviour. To describe 
fractionation it is necessary to define separate "daugh- 
ter" distributions p^'(cr) (i — 1,2...) which measure the 
distribution of the polydisperse attribute for each phase 
i. When the polydispersity is fixed, conservation of par- 
ticles implies that the weighted average of the daughter 
distributions equals the fixed overall density distribution, 
or "parent" p^(cr), For instance, at two-phase coexis- 
tence p(°)((r)=nf )/H = (1-£),o (1 V)+£p (2 V), with 
1— £ and £ the respective fractional volumes of the phases. 
This expression represents a generalisation of the Lever 
rule to polydisperse systems. 

To illustrate the central role of fractionation in polydis- 
perse phase behaviour, it is constructive to consider first 
the familiar case of the binodal curve of a monodisperse 
system in the density-temperature plane. This curve 
serves a dual purpose: on the one hand it describes the 
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range of overall densities for which phase coexistence oc- 
curs; and on the other hand it identifies the densities of 
the coexisting phases themselves. Now, for a polydis- 
perse system the range of overall (parent) densities that 
leads to coexistence is similarly delineated by a curve in 
n\°>-T space - the so-called "cloud" curve. However, the 
densities of the coexisting phases themselves do not in 
general coincide with the cloud curve. Instead, as one 
varies the parent density n' ' through the coexistence 
region at a fixed temperature (say), one generates an 
infinite sequence of pairs of differently fractionated co- 
existing phases. It is customary to single out the end 
points of this sequence for special attention, ie. the case 
of incipient phase separation that occurs when the value 
of nA°J coincides with the cloud curve. Under these con- 
ditions, one daughter phase has a fractional volume of es- 
sentially unity and consequently (from the Lever rule) a 
density distribution that is identical to the parent, while 
the other phase - known as the "shadow" - has an in- 
finitesimal fractional volume and a density distribution 
that deviates maximally from that of the parent. The 
curve formed by plotting the number density (or packing 
fraction) of the shadow phase as a function of tempera- 
ture is known as the shadow curve. 

Qualitative differences between the phase behaviour 
of monodisperse and polydisperse systems are also evi- 
dent in other projections of the full phase diagram. For 
instance, in the pressure-temperature plane, coexistence 
for a monodisperse system occurs along a simple line. 
By contrast for a polydisperse system, coexistence oc- 
curs within a region of the pressure-temperature phase 
diagram [H, Accordingly the task of exploring the 
full coexistence behaviour of a polydisperse system rep- 
resents a far greater endeavor than for a monodisperse 
system. 

Owing to the computational complexities associated 
with polydispersity, simulation studies of its effects on 
fluid-solid phase equilibria have been rather sparse. In 
pioneering work, Bolhuis and Kofke employed an isobaric 
semi-grand canonical ensemble approach to study the ef- 
fects of particle size polydispersity on hard sphere freez- 
ing [13, EH- Within this ensemble, Monte Carlo updates 
are performed that allow particles to change diameter, 
under the control of an imposed distribution of chemical 
potential differences p,(a\<7o) = p(a) — /x(cro), with <jq a 
reference particle diameter. The advantage of this ap- 
proach is two fold: it allows the instantaneous form of 
the density distribution p(cr) to fluctuate, thereby sam- 
pling many different realizations of the disorder in the 
course of a run, and it permits fractionation of the co- 
existing phases. However, in this early work no attempt 
was made to adapt the form of p,(cr) in order to ensure 
that the ensemble averaged density distribution p(o~) had 
a fixed functional form. Instead the activity distribu- 
tion exp[/3/2(<r)] was assigned a Gaussian form, peaked 
at ero, and various widths of the Gaussian were studied 
in order to change the degree of polydispersity. In such 
a fixed chemical potential representation, p(tr) can vary 



dramatically across the phase diagram. Additionally, co- 
existence occurs only along a line in the p — T plane (as 
in the monodisperse case) , rather than within a region as 
occurs experimentally. These artifacts limit the applica- 
bility of the results. 

To obtain the properties of the fluid-solid phase bound- 
ary, Bolhuis and Kofke employed Gibbs-Duhem integra- 
tion (GDI) techniques. This approach traces a coex- 
istence line in pressure-temperature space by integrat- 
ing its derivatives, obtainable from measurements of sin- 
gle phase properties near coexistence. The principal 
strengths of GDI is that it can quickly yield an estimate 
for a phase boundary and avoids mixed phases configura- 
tions. However, independent prior knowledge of a coex- 
istence state point is required in order to bootstrap the 
integration, and since there is no reconnection of the two 
phases beyond the initial starting point, integration er- 
rors can accumulate. The magnitude of these errors is 
difficult to quantify because there is no feedback mecha- 
nism to indicate when the integration has departed from 
the true coexistence line (provided one remains within 
the rather wide band of metastable states that borders 
this line). 

In more recent simulation work, Fernandez et al fl2| 
used a canonical ensemble approach to investigate the 
freezing of polydisperse soft spheres. However, owing to 
the small system size employed, the physical separation 
of the phases that would normally be expected to occur in 
a canonical ensemble simulation was (for the most part) 
absent and hence fractionation of the individual phases 
could not be measured. Consequently the "phase bound- 
ary" presented by Fernandez et al contains no informa- 
tion on the properties of any of the infinity of coexisting 
fluid-solid pairs, serving at best only as a rough estimate 
of the middle of the range of parent densities for which 
coexistence might occur. 

The method described in the present paper permits 
the accurate study of fluid-solid coexistence for systems 
of particles whose polydispersity is fixed. Accuracy is 
obtained by embedding the isobaric semi-grand canoni- 
cal ensemble within the PSMC framework, thus allowing 
the direct connection of the coexisting phases. Fixed 
polydispersity is realized by implementing an iterative 
reweighting scheme that adapts p,(a) and the pressure p 
such that the conjugate density distribution p(o~) matches 
a prescribed parent form, even within the coexistence 
region. Our paper is arranged as follows. In sec. |TT] 
we provide a summary of principal aspects of PSMC as 
applied to monodisperse systems, before proceeding to 
detail the extensions necessary to deal with polydisper- 
sity. Thereafter we describe how one can use the method 
to accurately determine fluid-solid coexistence proper- 
ties for systems having fixed polydispersity. Finally we 
provide some illustrative results for a system of polydis- 
perse soft spheres, obtaining cloud points, demonstrating 
the broadening of the coexistence region in the pressure- 
temperature plane, and quantifying fractionation effects. 
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II. PHASE SWITCH MONTE CARLO 

A full description of the statistical mechanical back- 
ground to PSMC as well as an in-depth discussion of 
its implementation for problems of fluid-solid coexistence 
has recently been given elsewhere (see ref. and shall 
not therefore be repeated here in full. Instead we shall 
focus somewhat more on the generalities of the method, 
as a prelude to describing the extensions required to deal 
with polydispersity. 

A. Monodisperse systems 

Let us work within an isobaric-isothermal (constant- 
NpT) ensemble, in which the particle number N , pressure 
p and temperature T are all prescribed [l3|. Then the 
relative stability of fluid (F) and crystalline solid (CS) 
phases is determined by the difference in their Gibbs free 
energies. This can in turn be obtained from the ratio of 
the a-priori probabilities of the phases [l|, Q : 

Ag ee g cs (N,p,T) - g F (N,p,T) = ^lnft F;C S • (1) 
with 



phase 7' simply by switching the set of reference sites 
{.R} 7 — » {-R} 7 , while holding the set of displacements 
{u} constant. This switch, which forms the heart of the 
method, can be incorporated in a global MC move (see 

fig.ru). 




Fluid Crystalline solid 



FIG. 1: Schematic illustration of the phase switch mecha- 
nism. The dots identify the representative sites {R} 1 in each 
phase; the displacement vectors {«} connect the centers of 
the distinguishable (numbered) particles to these sites. The 
switch operation shown swaps the representative sites of one 
phase for those of the other phase, whilst maintaining {u} 
constant. The particular configuration {u} shown is a "gate- 
way" state (see text) because the magnitude of the effective 
energy change under the switch is small. 



^F,CS 



P(F\N,p,T) 
P(CS\N,p,T) 



(2) 



In order to directly measure Tip cs > a MC procedure is 
required that samples both the solid and the fluid regions 
of configuration space in the course of a single simulation 
run at the prescribed values of N, p and T. From this 
one simply measures the probabilities P(F) and P(CS) 
of finding the system in each of the respective phases, 
and thence the probability ratio. Below we describe how 
PSMC facilitates such a measurement. 

PSMC takes as its starting point the specification of a 
reference configuration {R} 1 for each of the phases (la- 
beled 7) coexisting at the phase boundary. The specific 
choice of {i?} 7 is arbitrary, the only condition being that 
it should be a member of the set of pure phase config- 
urations identifiable as "belonging" to phase 7. For a 
crystalline phase, a suitably simple choice of {R} cs is 
the set of lattice sites; for a fluid, a suitable choice is a 
randomly chosen fluid configuration. 

The next step is to express the coordinates of each 
particle in phase 7 in terms of the displacement from its 
reference site, i.e. 



(3) 



Now, for displacement vectors that are sufficiently small 
in magnitude, one can clearly reversibly map any config- 
uration {r 7 } of phase 7 onto a configuration of another 



A complication arises however, because the displace- 
ments {u} typical for phase 7 will not, in general, be 
typical for phase 7'. Thus the switch operation will 
mainly propose high energy configurations of phase 7' 
which are unlikely to be accepted as a Metropolis up- 
date. This problem can be circumvented by employing 
extended sampling (biasing) techniques to seek out those 
displacements {u} for which the switch operation is en- 
ergetically favorable. These are the gateway configura- 
tions, which typically correspond to displacement vectors 
which are small in magnitude. Fig. [3 shows a schematic 
representation of the procedure. 

The requisite bias is administered with respect to an 
"order parameter" M. This macrovariable is defined such 
that the set of typical microstates (particle configura- 
tions) associated with its range, form a continuous phase 
space path linking the configurations of high statistical 
weight to the gateway configurations. In our formula- 
tion, the order parameter comes in two parts (or modes): 
'tether' and 'energy'. The tether mode serves to draw 
particles close to the representative sites to which they 
are nominally associated. Then, once all particles are 
within a prescribed distance of these sites, tether mode 
switches off and an energy mode order parameter takes 
over to guide particle to the gateway states for which 
the phase switch energy cost is sufficiently small to be 
accepted. 

To elaborate, let M mn denote the order parameter in 
mode to and phase 7. Then for tether mode we write 
to = T and define an associated order parameter 
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Fluid Crystalline solid 

FIG. 2: Schematic illustration of the phase switch operation 
in terms of the regions of configuration space associated with 
the fluid and crystalline solid phases. A bias (dashed ar- 
rows) is constructed such as to enhance the probability of 
the subsets of "gateway" states (the white islands) within the 
single-phase regions, from which the switch operation (the 
large arrow) will be accepted. Note that the switch accesses 
only one of the crystalline phase space replica fragments as- 
sociated with permutations of particles amongst lattice sites 



m t , 7 (M) = J2 max {°' "< - ' ( 4 ) 

where u,i = lu^/y 1 / 3 is the distance of particle i from 
its lattice site measured in units of the box length, and 
u c is a prescribed dimcnsionlcss threshold radius. Only 
particles whose displacement Ui exceeds this threshold 
contribute to Mt,t 

The tether mode is active iff vn > u c for at least one 
particle i, i.e. when Mr i7 > 0. Otherwise control hands 
over to the 'energy' mode m = £; its associated order 
parameter is defined by 

M £n {{u}) = sgn(A£ r7 ) ln(l + |A£y 7 |) , (5) 

where 

A£ Y7 = (E 1 - Ef) - (E Y - Ef) (6) 

measures the change (under the phase switch operation) 
of the enthalpy E y ({u},V) — <& 7 ({-u}) + pV with re- 
spect to its value £ in the representative microstate 
{i?} 7 , with the latter scaled to the instantaneous volume 
of phase 7 0, [l4[ . The presence of the logarithm in eq. \E\ 
is designed to moderate the scale of the contribution of 
the energy cost to Mg -y, which might otherwise become 
excessively large for particles with a strongly repulsive 
core to their interaction potential. 

Having defined a suitable order parameter, a biasing 
(weight) function is constructed to allow the system to 



reach the gateway states and hence - in the course of a 
sufficiently long run - switch back and forth repeatedly 
from fluid to solid. This weight function comes in four 
parts, one for each combination of phase and mode, and 
the interface between energy and tether parts are joined 
for continuity at their boundary. A number of techniques 
exist for constructing this weight function, but of those 
we have tested, we have found transition matrix Monte 
Carlo [HI to be (by far) the most efficient. Details of 
how to implement the transition matrix approach to ob- 
tain the weight function in the context of PSMC have 
previously been described elsewhere 0, H[ and we refer 
the interested reader to those articles for details. 

Turning now to the Monte Carlo procedure for sam- 
pling the order parameter distribution, in the context of 
monodisperse systems this involves four types of update, 
as outlined below. 

1. ' Particle translations'. A site (identified by one 
of the vectors in {R} cs or {R} ) is selected at 
random and a candidate state is chosen by incre- 
menting the position coordinate of the particle as- 
sociated with that site by a random vector whose 
components are drawn from a zero-mean uniform 
distribution of prescribed width. This operation 
changes both {r} and {u} 

2. 'Association swaps'. In this operation we choose 
two distinct sites at random (i and j say) and iden- 
tify the corresponding displacement vectors Ui and 
Hj. The candidate state is defined by the replace- 
ments 

Hi — ► Uj + Rj — Ri (7) 
Uj — > Ui + Ri — Rj (8) 

This process can be regarded as a change of as- 
sociation: the particle which was associated with 
site j is now associated with site i (and vice versa) . 
It changes the representation of the configuration 
(the coordinates {u}); but it leaves the physical 
configuration invariant. It need be applied only to 
the fluid phase where particles can wander far from 
their representative sites and need to be reined back 
in order to reach the gateway states. 

3. 'Volume moves'. The volume is also updated in 
the conventional way [13j, by a random walk of 
prescribed maximum step size, with particle posi- 
tion coordinates {r} and representative sites {i?} 7 
rescaled. 

4. 'Inter-phase switch'. The final type of MC update 
is the phase switch, which entails replacing one set 
of the representative configuration vectors, {i?} 7 
say, by the other, {-R} 7 . The switch should also 
incorporate an appropriate volume scaling of the 
system. 
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Details of the appropriate biased acceptance probabil- 
ities for these moves can be found in the original articles 

Using these MC updates, together with a suitable esti- 
mate of the weight function, the sampling process can be 
initiated. During this process one accumulates (initially 
in the form of a list [l||) the joint probability distribution 
of the order parameter, the volume, the energy and the 
phase label 7. Subsequently the effects of the imposed 
bias are unfolded from this distribution (in the standard 
manner flSH) to provide a direct measure of the rela- 
tive probabilities of the two phases: 

= JdM dVd£ P(M,V,£,F) 
F ' cs / dM dV d£ P(M, V, £, CS) ' 1 ' 

from which the Gibbs free-energy-density difference fol- 
lows directly via Eqs.[T]and[21 In practice these integrals 
are estimated by simple sums over the list of sampled val- 
ues of M, V, £ , 7. Use of histogram reweighting (HR) [l7| 
permits exploration of values of the p and T in the neigh- 
bourhood of a given simulation state point and serves as 
an invaluable aid for pinpointing the coexistence param- 
eters at which the free energy difference vanishes. 

B. Extension to polydisperse systems 

The natural choice of simulation ensemble for the study 
of F-CS coexistence in polydisperse s yste ms is the iso- 
baric, semi-grand-canonical ensemble [181 ]. Within this 
ensemble, the particle number N, pressure p, temper- 
ature T, and a distribution of chemical potential differ- 
ences p(o~\o~o) are all prescribed, while the system volume 
V, the energy, and the form of the instantaneous density 
distribution p(cr) all fluctuate. The fluctuation in p(o~), 
although subject to the constraint V J p(a)da = N, nev- 
ertheless permits the sampling of many realizations of 
the polydisperse disorder, thus ameliorating finite-size 
effects. Moreover, in conjunction with volume fluctua- 
tions, it facilitates separation into differently fractionated 
phases. 

Operationally, the sole difference between the isobaric, 
semi-grand-canonical ensemble and the constant- A^T en- 
semble used for PSMC is that one implements MC up- 
dates that select a particle at random and attempt to 
change its diameter a by a random amount drawn from 
a zero-mean uniform distribution. This proposal is ac- 
cepted or rejected with a Metropolis probability con- 
trolled by the change in the internal energy and chemical 
potential [Hj]: 

p acc = min [1, exp (-/3[A$ + ~p(o) - £((/)])] , 

where A<& is the internal energy change associated with 
the resizing operation. Accordingly, one can readily ex- 
tend the PSMC framework to deal with polydispersity by 



specifying a form for //(er|ero) and supplementing the four 
Monte Carlo operations described in Sec. Ill Al with a re- 
sizing move. For the purposes of histogram reweighting, 
one samples the joint probability distribution of the order 
parameter, the volume, the energy and the instantaneous 
density distribution p(o~). Again this can be conveniently 
stored as a simple list of sampled values. 

III. FIXED POLYDISPERSITY 

For simulations of fixed polydispersity at some given 
N and T, one requires that both the pressure p and p(a) 
are such that a suitably defined ensemble-averaged den- 
sity distribution matches the prescribed parent p^ (a) = 
n^°'f(a). Unfortunately, the task of determining the req- 
uisite p and p(a) is severely complicated by the fact they 
are an unknown functional of the parent. Essentially, 
therefore, one is faced with solving an inverse problem 

To deal with this problem we employ a scheme similar 
to one recently proposed in the context of grand canonical 
ensemble studies of polydisperse phase coexistence [20|. 
The method relies on the fact that data accumulated at 
a given set of p and p(a) can be used to extrapolate 
(via histogram reweighting) to other nearby sets of p and 
/2(er), without further simulation. In our description we 
shall assume for convenience that while the particle size 
o~ is to be treated as a continuous variable, distributions 
defined on a, such as p(cr), and /i(cr), are represented as 
histograms formed by discretising the range of a into a 
finite number of bins. 

One proceeds as follows. The first step is to decide 
which of the infinity of coexisting pairs of phases inside 
the cloud curve one would like to determine for a given 
temperature T . To do so one specifies a value for the 
parameter £ which measures the fractional volume of the 
solid phase, and thus parameterizes the dilution lines be- 
tween the liquid phase cloud (£ = 0) and the solid phase 
cloud (£ = 1). The general strategy for determining the 
coexistence properties at the given £ (and T) is then to 
tune p, p(o~) and iteratively, such as to simultane- 
ously satisfy both a generalized lever rule and equality of 
the a-priori probabilities of the phases, ie. 

= (1-0P (1 V)+£P (2 V), (10a) 
^f,CS = 1- (10b) 

In the first of these constraints (Eq. UOap , the ensem- 
ble averaged daughter density distributions p^(a) and 
p^ (a) are assigned by averaging only over configurations 
belonging to the fluid and solid phases respectively. The 
deviation of the weighted sum of the daughter distribu- 
tions p(a) = (1 — 0^ 1 H £T ) + iP < ~ 2 H (J ) lrom the target 
n (o) j(o)(ct) is conveniently quantified by a 'cost' value: 

A = J \p(a)-nWfW(a)\da. (11) 
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Requiring 7?.p,CS = 1 in the second constraint (Eq. I10b| . 
ensures that finite-size errors in coexistence parameters 
are exponentially small in the system volume [2(| HH • 

The iterative determination of p, jl(cr) and such as 
to satisfy Eqs. UOal and 1 1 Obi then proceeds thus: 

1. Guess an initial value of the parent density n^ ' cor- 
responding to the chosen value of £. (If one starts 
with a small degree of polydispersity, a suitable es- 
timate can be obtained from the densities of the 
coexisting phases in the known monodisperse limit, 
together with the lever rule.) 

2. Tune the pressure p (within the HR scheme) such 
as to minimize A. 

3. Similarly tune p(<r) (within the HR scheme, see be- 
low) such as to minimize A. 

4. Repeat steps 2 and 3 until A < Tolerance. 

5. Measure the corresponding value of 7?-f,CS- 

6. if |7?-f ; cs — 1| < Tolerance, finish, otherwise vary 

and repeat from step 2. 

The minimization of A with respect to variations in 
p (step 2) can be easily automated using standard 1- 
dimensional minimization algorithms such as the "brent" 
routine described in Numerical Recipes [22j. The same 
applies to the minimization of {IZf^cs ~ 1| with respect 
to variations in in step 6. In step 3 the minimization 
of A with respect to variations in jj,(cr) is most readily 
achieved [23| using the following simple iterative scheme 
for p,(a): 



/2 fc+ i(cr) = jSfc(o-) +a\n y — J , (12) 

for iteration k — ► k + 1. This update is applied simul- 
taneously to all entries in the histogram of /i(c), and 
thereafter the distribution is shifted so that /2(er ) = 0. 
The quantity < a < 1 appearing in eq. [12] is a damping 
factor, the value of which may be tuned to optimize the 
rate of convergence. 

In practice we find that the minimization of A with 
respect to variations in p{a) operates well provided p(a) 
is sufficiently close to f(°\a) for histogram reweight- 
ing to be effective. Note that (as described in [2(|) it is 
important that in steps 4 and 6, one iterates to a very 
high tolerance in order to ensure that the finite-size ef- 
fects are exponentially small in the system size. Typically 
we iterated until both A and \TZp cs — 1| were less than 
1(T 12 . 

The values of nS ' and p resulting from the application 
of the above procedure are the desired parent density and 
pressure corresponding to the nominated value of £. For 
£ = and £ = 1, these are (respectively) the liquid and 
solid cloud point densities and pressures, while shadow 



points are given by the properties of the coexisting incip- 
ient daughter phase, which may be simply read off from 
the appropriate peak positions of the cloud point dis- 
tributions of quantities such as the fluctuating packing 
fraction. 

Finally in this section, we note that whilst in the in- 
terests of clarity we have not included a description of 
temperature reweighting in our procedure, it is straight- 
forward to incorporate such. Doing so allows one to scan 
the p—T plane in a stepwise fashion, and thus, ultimately, 
determine the entire transition region. 

IV. SOLID-FLUID COEXISTENCE OF 
SIZE-DISPERSE SOFT SPHERES 

We have applied our methodology to obtain the cloud 
and shadow point properties of a polydisperse system of 
softly repulsive spheres described by the r~ 12 potential 

v(r)=e(a ij /r) 12 , (13) 

with Uij = (<Tj + Oj)/2. In fact, for a given choice of /(cr) 
the thermodynamic state of this model does not depend 
on n^ ) and T separately but only on the combination 
n^^e/fcsT) 1 / 4 . Thus, the coexistence points for various 
temperatures scale exactly onto one another [24[ and one 
therefore only needs to determine coexistence properties 
for a single temperature. Accordingly in this work we 
shall consider the case (e/ksT) = 1. 

In our simulations, the interparticle potential was trun- 
cated at half the box size and periodic boundary condi- 
tions were applied. Since for the system sizes studied the 
value of the potential is extremely small at the typical 
cutoff radius, no correction was applied for the trunca- 
tion. 

We chose to study a parent distribution of the top-hat 
form defined by 

, f (2c)- 1 ifl-c<a<l + c (u] 
^ ' \ otherwise ' ' 

where (without loss of generality), the mean particle di- 
ameter has been set to a — 1.0. In this initial study 
we report results for a system of TV = 256 particles and 
rather narrow parent distributions c < 0.05, for which 
the dimensionless degree of polydispersity (ratio of stan- 
dard deviation to mean of /(ct) is S = c/y/3 < 3%. For 
the purposes of forming p(a) and pier), each unit of a 
was discretized into 500 bins. 

As a preliminary step we determined the coexistence 
parameters of the transition from fluid to face-centred- 
cubic (fee) solid in the monodisperse limit. This was 
done using the standard formulation of PSMC (see 
sec. IIIAI) with the results [6]: p CO ex — 22.32(3), pf = 
1.148(9), pes = 1.190(9). Thereafter we attempted to 
locate the fluid phase cloud point (£ = 0) for a narrow 



7 



top-hat parent having c = 0.01. To this end we initialized 
the chemical potential difference distribution as p.(cr) = 
(for 0.99 < a < 1.01) and fj,(a) = -100 otherwise, and 
assigned the pressure the value p — 22.32 pertaining to 
the monodisperse limit. We then performed a long PSMC 
run, the results of which were reweighted (using the pro- 
cedure described in Sec. IIIII together with the monodis- 
perse value n^ ' = 1.148 as the initial guess for the fluid 
cloud point density), to yield accurate estimates of the 
fluid phase cloud point pressure, parent density and 
chemical potential difference distribution p.(<r), as well as 
the shadow phase daughter distribution. 

In order to progress to higher degrees of polydispersity, 
we proceeded in a stepwise fashion. The form of jl(a) pre- 
viously determined for c = 0.01 was extrapolated via a 
quadratic fit to cover the range 0.98 < a < 1.02, while 
the pressure for c = 0.02 was estimated by linearly ex- 
trapolating the results for c = 0.00 and c = 0.01. A new 
PSMC run was then performed, the results of which were 
reweighted to give accurate estimates of the cloud point 
parameters for the top hat parent having c = 0.02. In 
this manner we were able to steadily increase the degree 
of polydispersity, measuring the cloud point pressure and 
parent density as we went. In an identical fashion we de- 
termined the dependence on polydispersity of the solid 
cloud point parameters, by setting £ = 1 in the above 
procedure. 




0.02 
0.05 



PF 



,(<>) 



PCS 



22.32(3) 1.148(9) 22.32(3) 1.190(9) 
22.46(3) 1.149(1) 22.49(3) 1.192(2) 
23.21(4) 1.159(1) 23.39(5) 1.201(2) 



TABLE I: Estimates of the coexistence pressure and parent 
density at the fluid and fcc-solid cloud points for the half- 
widths c of the top-hat parent distribution indicated. 

Table. Q] shows the fluid and solid cloud point param- 
eters for c = 0,0.02 and 0.05. We note that our results 
for the monodisperse limit (c = 0) are consistent with 
existing estimates for the same model obtained via ther- 
modynamic integration [24j. Our data for c > exhibit 
a clear separation of the cloud point pressures as poly- 
dispersity increases, reflecting the polydispersity-induced 
broadening of the coexistence region described in Sec. [U 
Also apparent is an increase in the cloud point densities 
with increasing polydispersity, a finding which is in ac- 
cord with theoretical predictions for polydisperse hard 
spheres [25]. In Fig. [3] we plot the form of the shadow 
phase daughter distributions for c — 0.05. This shows 
that at the fluid phase cloud point the incipient solid 
phase contains a majority of large particles compared to 
the parent, whilst at the solid phase cloud point, the 
incipient (liquid) phase contains a majority of small par- 
ticles. This too is in accord with the predictions for mod- 
els of polydisperse hard spheres [25[ . Fig. [3^b) shows the 
corresponding cloud point forms of the chemical potential 
difference distribution /2(<r). 
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□ Fluid shadow 

— Fluid cloud (5=0) 

— Solid cloud (|=1) 



1.94 
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FIG. 3: (a) Estimates of the cloud and shadow point density 
distributions for c == 0.05. (b) Corresponding distribution 
of chemical potential differences, for <7o = 0.95. Statistical 
uncertainties do not exceed the symbol sizes. 



Finally, we consider how one obtains the properties 
of shadow phases such as their average packing fraction 
or average number density. This is illustrated in Fig. 0] 
which shows (for c = 0.05) the form of the probability dis- 
tribution of the fluctuating packing fraction at the solid 
cloud point (£ = 1), ie. the distributions P(r)^') and 
P(r)W) with r?W = (tt/6) J a 3 p^(a). The peak at large 
rj corresponds to the solid cloud phase and the peak at 
lower rj corresponds to the fluid shadow phase. The aver- 
age of the shadow point packing fraction fj can be simply 
read off from the position of the shadow phase peak. For 
a system with a non-trivial temperature dependence of 
the phase boundary, such plots (and analogous one for 
the number density), obtained for a range of tempera- 
tures, allow construction of the shadow curve in cither 
the packing fraction or number density representation. 



V. CONCLUSIONS 

In summary, we have presented an accurate approach 
for determining the fluid-solid coexistence properties of 
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0.61 0.62 0.63 

volume fraction r\ 



0.65 



FIG. 4: The distribution of the fluctuating packing fraction 
P(rf) in the solid cloud phase and in the shadow daughter 
phase for c = 0.05, as described in the text. Statistical uncer- 
tainties are comparable with the symbol sizes. 



polydisperse fluids under the experimentally relevant 
constraint that the parent distribution has a fixed func- 
tional form. The method is capable of determining the 
properties of any of the infinity of coexistence state points 
that pertain inside the cloud curve. By illustration we 
have determined the cloud and shadow properties at the 
freezing transition for a system of soft spheres having 
small degrees of polydispersity. 



Like the bare PSMC method on which it is based, our 
approach requires a significant investment of computa- 
tional effort in order to reap the benefits of accurate es- 
timates of coexistence properties. The results presented 
here required about 4 weeks of CPU time on an 8-core, 
3GHz processor. Nevertheless, it should be straightfor- 
ward to extend our study to higher degrees of polydisper- 
sity and larger system sizes, albeit with a correspondingly 
increased computational expenditure. 

Thus our method could prove useful in helping to inves- 
tigate interesting theoretical predictions [2 51 ] concerning 
the fluid-solid coexistence properties of polydisperse flu- 
ids. Specifically, it has been proposed that while on the 
fluid cloud curve the degree of polydispersity of the fluid 
parent can be increased to rather large values of 8 > 14%, 
in the coexisting (shadow) solid daughter phase, 8 never 
exceeds the more modest value 8 w 6%. On the other 
hand, on the solid phase cloud curve, the single solid 
parent is itself predicted to become unstable at about 
8 w 6% and phase separates at a triple point into two 
solid phases and a liquid. We hope to be able to address 
these issues in future work. 
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